function [corrplanes, sups, corratts, corrpoints]=CARE(data, maxs, fmax,  minsup)
[M,N] = size(data);
candiatts = [1:N];
corrplanes={};
corratts={};
corrpoints={};
sups=[];

subatts = [1,2,3];
[subcorrplanes, subsups, subcorratts, subcorrpoints]=findlocalcorrelations(data, maxs, subatts, fmax, minsup);
corrplanes ={corrplanes{1:end},subcorrplanes{1:end}};
corratts ={corratts{1:end},subcorratts{1:end}};
corrpoints={corrpoints{1:end},subcorrpoints{1:end}};
sups=[sups, subsups];
    
% for i=1:N
%     subatts = [i];
%     [subcorrplanes, subsups, subcorratts]=findlocalcorrelations(data, maxs, subatts, fmax, minsup);
%     corrplanes ={corrplanes{1:end},subcorrplanes{1:end}};
%     corratts ={corratts{1:end},subcorratts{1:end}};
%     sups=[sups, subsups];
% end

function [corrplanes, sups, corratts,corrpoints]=findlocalcorrelations(data, maxs, attrs, fmax, minsup)
[M,N] = size(data);
corrplanes={};
corratts={};
corrpoints={};
subdata = [];
sups =[];
if length(attrs) > maxs
    return;
end
for i=1:length(attrs)
    subdata = [subdata, data(:,attrs(i))];
end
disp(attrs);
[plane, numpoint, f, points]=primaryplane(subdata, fmax, minsup);
sup = numpoint/M;
if f >= -eps && f <= fmax && sup >= minsup
    corrplanes{end+1}=plane;
    corratts{end+1}=attrs;
    corrpoints{end+1}=points;
    sups = [sups, sup];
end
if attrs(end)+1 > N
    return;
end
candiatts = (attrs(end)+1):N;
for i=1:length(candiatts)
    subatts=[attrs, candiatts(i)];
    [subcorrplanes, subsups, subcorratts]=findlocalcorrelations(data, maxs, subatts, fmax,minsup);
    corrplanes ={corrplanes{1:end},subcorrplanes{1:end}};
    corratts ={corratts{1:end},subcorratts{1:end}}; 
    corrpoints={corrpoints{1:end},subcorrpoints{1:end}};
    sups = [sups, subsups];
end

% Find the primary plane
function [plane, numpoint, f, points]=primaryplane(data, fmax, minsup)
tmpdata = data;
[M,N]=size(data);
[M2,N]=size(tmpdata);
points=[];
while M2>=M*minsup
    [pri,f]=smallestPC(tmpdata);
    cmean=mean(tmpdata);
    const = -pri'*cmean';
    if f<=fmax
        plane=[pri',const];
        numpoint=M2;
        points=tmpdata;
        return;
    end
    maxDist=0;
    maxDistPI=-1;
    for i=1:M2
        dist = abs(pri'*tmpdata(i,:)'+const)/norm(pri');
        if dist>maxDist
            maxDist = dist;
            maxDistPI = i;
        end
    end
    M2=M2-1;
    tmpdata=[tmpdata(1:maxDistPI-1,:);tmpdata(maxDistPI+1:end,:)];
end
plane=[];
numpoint=0;
f=0;

% Find the smallest component vector
function [pri,f]=smallestPC(data)
[M,N]=size(data);
if N == 1
    pri = [1.0];
    f=1.0;
    return
end
mn=mean(data,1);
newdata=data - repmat(mn,M,1);
covariance=1/(M-1)*newdata'*newdata;
[PC,V]=eig(covariance);
V=diag(V);
[junk, rindices]=sort(V);
pri=PC(:,rindices(1));
if sum(V) == 0
    f = 0;
else
    f=V(rindices(1))/sum(V);
end

